;
; Model formula tools for XLISP-Stat regression methods
; Modified and adapted from T. Lumley
; Pedro Valero 2003
;Note arguments for excluded categories are unused yet
(defun design-matrix (matrix types model vars &optional excluded-categories)
  (let* (
         (data-matrix-column-list (column-list matrix))
         (variables-list vars)
         (types types)
         (res (mapcar #'(lambda (el)
                          (cond 
                            ((listp el)
                             (let ((vars (mapcar #'(lambda (v) (coerce v 'list))
                                                 (select data-matrix-column-list el)))
                                   (is-factor (mapcar #'(lambda (l) 
                                                          (equal "CATEGORY" 
                                                                 (string-upcase l)))
                                                      (select types el)))
                                   (namebases (select variables-list el))
                                   )
                               
                               (interaction vars
                                            :is-factor is-factor
                                        :namebases namebases)
                               ))
                        ((equal "CATEGORY" (string-upcase (select types el)))
                         (let ((var (coerce (select data-matrix-column-list el) 'list))
                               )
                           (factor  var
                                    :namebase (select variables-list el))
                           ))
                            ((equal "NUMERIC" (string-upcase (select types el)))
                         (let ((var (coerce (select data-matrix-column-list el) 'list))
                               (variable-name (select variables-list el))
                               )
                           (term var
                                 :namebase variable-name))
                             )))
                  model))
         )
    (setf blocks (mapcar #'(lambda (b) (if (matrixp b) (array-dimension b '1) 1)) (mapcar 'first res)))
    (list (apply 'bind-columns (mapcar 'first res)) (combine (mapcar 'second res)) blocks)))


(defun factor (x &key namebase (excluded-category nil)) "Defines treatment contrast matrix for x and optionally a set of names based on namebase"
  (let* (
	 (xlist (remove-duplicates (coerce x 'list) :test #'equalp))
         (xlist (if excluded-category 
                    (combine excluded-category (sort-data (remove excluded-category xlist)))
                    (sort-data xlist)))
       	 (p (- (length xlist) 1))
	 (rows (row-list (select (identity-matrix (+ p 1)) (iseq 0 p) (iseq 1 p))))
	 (decoder (mapcar #'cons xlist rows))
	 (xmatrix (apply #'bind-rows (mapcar #'(lambda (xx) (cdr (assoc xx decoder :test #'equalp))) x)))
	 (xnames (if (null  namebase) nil (cdr (mapcar #'(lambda (xx) (princ-to-string (list namebase xx) )) xlist))))
    )
    (if (null namebase) xmatrix (list xmatrix xnames))
)
)


(defun interaction ( xs &key (is-factor (repeat t (length xs))) namebases ) "Interactions of any order for categorical and continuous variables"
  (let*  ( (n (length (first xs)))
	  (zlist nil)
	  (znames nil)
	  (p (length xs))
	  (names (if (null namebases) (repeat "." p) namebases))
	  (junk (dotimes (i p)
			(let* (
			       (thisz  (if (elt is-factor i) (factor (elt xs i) :namebase (elt names i)) (list (bind-columns (elt xs i)) (list (elt names i)))))
			       )
			  (if (null zlist) (setf zlist (column-list (first thisz)))   (setf zlist (apply #'append (mapcar #'(lambda (xx) (mapcar #'(lambda (yy) (*  xx yy)) (column-list (first thisz)))) zlist))))
			  (if (null znames) (setf znames (second thisz))  (setf znames (apply #'append (mapcar #'(lambda (xx) (mapcar #'(lambda (yy) (princ-to-string (list xx yy) )) (second thisz))) znames))))
			  )
			)
		)
	  (zmat (apply #'bind-columns zlist))
	  )
    (if (null namebases) zmat (list zmat znames))
))
	  
			  

;
(defun block-test (index beta covmat &key blockname names (block-only nil))
  (let* (
	 (subbeta (select beta index))
	 (subcov  (select covmat index index))
	 (waldchisq (inner-product subbeta (matmult (inverse subcov) subbeta)))
	 (waldp (- 1 (chisq-cdf waldchisq (length index))))
	 (subse (sqrt (diagonal subcov)))
	 (subz (/ subbeta subse))
	 (subp (* 2 (- 1 (normal-cdf (abs subz)))))
	 (nn (if (null names) (repeat "" (length index)) (select names index)))
	 (blockn (if (null blockname) "block" blockname))
	 )
    (format t "~a~20t~13,5g~35t~,4f~%" blockn waldchisq waldp)
    (cond ((null block-only)
	(format t "~5t Variable~25t Estimate~40t Std.Err.~55t p-value~%")
	(dolist (i (iseq 0 (- (length index) 1)))
	(format t "~5t~a~25t~13,5g~40t(~,6g)~55t~,4f~%" (select nn i) (select subbeta i) (select subse i) (select subp i)))
	)
      )
    )
)


(defun design (varlist) (apply #'bind-columns (mapcar #'first varlist)))

(defun names (varlist) (apply #'append (mapcar #'second varlist)))

(defun term (x &key namebase) "A metric predictor variable"  (list x namebase))

(defmeth scatterplot-proto :my-new-plot (x &optional (y nil) )
  (let*
    ((num-points (send self :num-points))
     (color (send self :point-color (iseq num-points)))
     (symbol (send self :point-symbol (iseq num-points)))
     (state (send self :point-state (iseq num-points)))
     (selection (send self :selection))
     (showing (send self :point-showing (iseq num-points)))
     (labels (send self :point-label (iseq num-points)))
     )
    (send self :start-buffering)
    (send self :clear-points)
    
    (if y (send self :add-points x y)
        (send self :add-points x))
    (send self :point-color (iseq num-points) color)
    (send self :point-symbol (iseq num-points) symbol)
    (send self :point-state (iseq num-points) state)
    (send self :selection selection)
    (send self :point-label (iseq num-points) labels)
    (send self :redraw)
    (send self :buffer-to-screen)
     ))